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ABSTRACT 


We present the first measurement of the mass function of free-floating planets (FFP) or very wide orbit planets down 
to an Earth mass, based on microlensing data from the MOA-II survey in 2006-2014. The shortest duration event has 
an Einstein radius crossing time of tg = 0.057 + 0.016 days and an angular Einstein radius of 0g = 0.90 + 0.14 yas. 
There are seven short events with tg « 0.5 day, which are likely to be due to planets. The detection efficiency for short 
events depends on both £g and Og, and we measure this with image-level simulations for the first time. These short 
events can be well modeled by a power-law mass function, dN4/dlog M = (2.18*9 22) x (M/8 Ma)-^* dex~'star~! 
with o4 = 0.96*057 for M/Mọo < 0.02. This implies a total of f = 21*75 FFP or wide orbit planets in the mass 
range 0.33 < M/Ma < 6660 per star, with a total FFP mass of m = 80* 13 Ma per star. The number of FFP is 1g 
times the number of planets in wide orbits (beyond the snow line), while the total masses are of the same order. This 
suggests that the FFPs have been ejected from bound planetary systems that may have had an initial mass function 
with a power-law index of o ~ 0.9, which would imply a total number of novus planets star~! and a total mass of 
ITL 5 Mọ star™!. 


Keywords: gravitational microlensing; exoplanet; Free floating planets 
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1. INTRODUCTION 


Gravitational microlensing observations toward the 
Galactic bulge (Galactic Bulge) enable exoplanet 
searches (Mao & Paczyński 1991; Gaudi et al. 2008; 
Bennett et al. 2010; Suzuki et al. 2016; Koshimoto et 
al. 2021b), and the measurement of the stellar and sub- 
stellar mass functions (MFs) (Paczyński 1991; Sumi et 
al. 2011; Mróz et al. 2017, 2019, 2020a). 

Sumi et al. (2011) first interpreted the detection of 
short timescale (0.5 < tg/day < 2) microlensing events 
as evidence for the existence of a population of free- 
floating planets (FFP) and/or wide orbit planets. While 
that analysis was limited by the small number of events 
found in a 2 year subset of the survey by the Microlens- 
ing Observation in Astrophysics (MOA) group (Sumi 
et al. 2003) in collaboration with Optical Gravitational 
Lensing Experiment (OGLE) (Udalski et al. 1994), it 
opened up the field of FFP studies using microlensing. 

Mróz et al. (2017) extended the work by using a larger 
sample from 5 years of the OGLE survey. They discov- 
ered 6 events with timescales shorter (tg ~ 0.2 day) than 
those in the previous work. These events are separated 
from the longer events by a gap around tg ~ 0.5 day 
which implying the possibility of a several Earth-mass 
FFP population. 

These studies are based on distribution of tg, in which 
tg is proportional to the square root of the lens mass M 
as follows, 
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Here, k = 4G/(c?au) = 8.144mas/Mg and we expect 
tg ~ 0.1 day assuming typical value of the lens-source 
relative parallax: mq = 7; /—7; | —lau(Dj!-D;!)- 
18yas for the bulge lens and a typical value of the lens- 
source proper motion in the direction of the Galactic 
center of fre) = 5 masyr™!. The lens mass M, the 
distance D to the lens and the relative proper motion 
Ure] are degenerate in the observable tg. This means 
that the mass function of the lens population has to be 
determined statistically, assuming a model of the star 
population density and velocities in the Galaxy. 

Mróz et al. (2018) found the first short (tg — 0.32 
day) event showing the Finite Source (FS) effect, i.e., a 
finite source and a single point lens (FSPL), in which 
one can measure a FS parameter p = 0./0g. Here 6, is 
the angler source radius which can be estimated from an 
empirical relation with the source magnitude and color. 


tg 


The 0g is an angular Einstein radius given by 


On = zam VKM Trel- (2) 


tg 


This value of 0g can give us an inferred mass of the 
lens with better accuracy as we can eliminate one of the 
three-fold degenerate terms which affecte tg, namely, 


Hrel: 
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So far, seven short FSPL events have been discov- 
ered (Mróz et al. 2018, 2019b, 2020b,c; Kim et al. 2021; 
Ryu et al. 2021). All of these have 0g < 10 yas, imply- 
ing that their lenses are most likely of planetary mass. 
All of these sources are red giants because their angu- 
lar radii, i.e., cross-section, are significantly larger than 
main sequence stars. 

Mróz et al. (2020b) found the short FSPL event, 
OGLE-2016-BLG-1928, with the smallest value of 0g = 
0.842 + 0.064uas to date. Its lens is the first terrestrial 
mass FFP candidate and the first evidence of such pop- 
ulation. 

Kim et al. (2021) began a new approach to probing 
the FFP population by focusing on analyzing the 0g 
distribution in events with giant sources. Ryu et al. 
(2021) found a gap at 10 < 0g/uas < 30 in the cumu- 
lative 0g distribution, which suggests a separation be- 
tween the planetary mass population and other known 
populations, like brown dwarfs. 

Gould et al. (2022) completed the analysis of 29 
FSPL giant-source events found in the 2016-2019 KMT- 
Net survey. They presented the Og distributions down 
to 0g = 4.35 uas and confirmed that there is a clear gap 
in the distribution of 0g at 9 < 0g/uas < 26. They 
note that it is consistent with the gap in the tg distri- 
bution shown by Mróz et al. (2017), indicating the ex- 
istence of the low mass FFP population. They modeled 
the 0g distribution with a power law MF for the FFP 
and found dNppp/dlog M = (0.4 + 0.2)(M/38Mg) ? 
dex-!star-!, with 0.9 « p « 1.2. This implies that the 
number of FFPs is at least an order of magnitude larger 
than that for known bound planets. 

In this paper, we present the distributions 0g and tg 
values for the microlensing events toward the Galac- 
tic Bulge from 9 years of the MOA-II survey. We also 
present the MF of the planetary mass objects using the 
tg distribution. We describe the data in section § 2. We 
show the 0g distribution in section 83. We present the 
tg distribution and the best-fit MF in 84. The discus- 
sion and conclusions are given in section 85. 
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2. DATA 


We use the microlensing sample selected from the 
MOA-II high cadence photometric survey toward the 
Galactic Bulge in the 2006-2014 seasons (Koshimoto et 
al. 2023, hereafter K23). MOA-II uses the 1.8-m MOA- 
II telescope which has a 2.18 deg? field of view (FOV) 
and which is located at the Mt. John University Obser- 
vatory, New Zealand!. 

K23 used an analysis method similar to what was used 
by Sumi et al. (2011, 2013), but includes a correction 
of systematic errors and takes into account the finite 
source effect. They applied a de-trending code to all 
light curves to remove the systematic errors that corre- 
late with seeing and airmass due to differential refrac- 
tion, differential extinction and relative proper motion 
of stars in the same way as in Bennett et al. (2012) and 
Sumi et al. (2016a). These corrections are important as 
they result in higher confidence in the light curve fitting 
parameters. 

K23 selected light curves with a single instantaneous 
brightening episode and a flat constant baseline, which 
can be well fit with a point-source point-lens (PSPL) 
microlensing model (Paczyński 1986). In addition to 
PSPL, they modeled the events with a FSPL model 
(Bozza et al. 2018), which is especially important for 
short events. These are the major improvements com- 
pared to the previous analysis in Sumi et al. (2011, 2013) 
in addition to the extension of the survey duration. 

Although they identified 6,111 microlensing candi- 
dates, they selected only 3,554 and 3,535 objects as the 
statistical sample using the two relatively strict criteria 
CR1 and CR2, respectively. Here, CR2 was defined as 
the stricter criteria compared to their nominal criteria 
CR1 to check the effect of the choice of the criteria on a 
statistical study. These strict criteria ensure that tg is 
well constrained for each event and reject any contami- 
nation. 

Sumi et al. (2011) reported 10 short events with tg « 2 
days in the 2006-2007 dataset. Only 5 and 4 events sur- 
vived following the application of CR1 and CR2, respec- 
tively. This is because the fitting results changed due to 
the re-reduction of the dataset. On the other hand, two 
events are newly found resulting 7 and 6 events follow- 
ing the application of CR1 and CR2, respectively. As a 
result, the excess at tg = 0.5 — 2 day in the tp distribu- 
tion is not significant anymore, however, an even shorter 
event MOA-9y-6057 (tg = 0.22 + 0.06 day) is added. 


3. ANGULAR EINSTEIN RADIUS DISTRIBUTION 


l nttps://www.massey.ac.nz/-iabond/moa/alerts/ 
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Figure 1. Observed cumulative distribution of 0g for 13 
FSPL events from MOA (red line) and 29 FSPL events from 
KMTNet (black line) (Gould et al. 2022). The blue line in- 
dicates 0g = 0.8424- 0.064 of terrestrial mass FFP candidate, 
OGLE-2016-BLG-1928 (Mróz et al. 2020b). 


There are 13 FSPL events with 0g measurements in 
the sample, including two FFP candidates, MOA-9y- 
5919 and MOA-9y-770, that have terrestrial and Nep- 
tune masses, respectively. See K23 for the light curves 
and detailed parameters of the 13 events. 

The red line in Figure 1 indicates the cumulative dis- 
tribution of Og from Table 7 of K23. The black line 
indicates the distribution of 29 FFPs by Gould et al. 

(2022) normalized to 13 events as a comparison. Al- 
though these can not be directly compared because these 
are not corrected for detection efficiencies, the general 
trends seen Figure 1 may give us some insights. 

The distributions are consistent for 0g > 30 pas, 
where the effect of the detection efficiencies are likely 
small. There is a gap around 5 < 0g/pas < 70 which 
is roughly consistent with the gap at 10 < 0g/puas < 30 
found by Ryu et al. (2021) and Gould et al. (2022). 
'This gap confirmed the existence of the planetary mass 
population as distinct and separated from the stel- 
lar/brown dwarf population as indicated by Gould et 
al. (2022). 

'The MOA cumulative distribution shows fewer events 
over 30 < ÓOg/uas < 70 compared to Gould et al. 
(2022). This may be just due to the small number of 
statistics. But note that K23 found a brown dwarf can- 
didate MOA-9y-1944 with 0g = 46.1 10.5 uas although 
this is not in the final sample for statistical analysis be- 
cause the source magnitude of Is = 21.91 mag is fainter 
than the threshold of Is < 21.4 mag. 

In our sample, there is one event with a very small 
value of 0g of 0.90+0.14 pas. This confirms the existence 
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Table 1. Comparison of parameters of short FS events with known FFP candidates. 


field-chip-sub-ID tg p Is. 0.. On reference 
(day) (mag) (pas) (uas) 
MOA-9y-5919 0.057 + 0.016 1.40 + 0.46 17.23 1.26 +0.48 0.90 + 0.14 K23 
MOA-9y-770 0.315 + 0.017 1.08 + 0.07 14.71 5.13 +0.86 4.73 + 0.75 K23 
OGLE-2016-BLG-1928 0.0288 *00024 3.397079 15.78 2.85 + 0.20 0.842 0.064 Mróz et al. (2020b) 
KMT-2019-BLG-2073 0.272 + 0.007 1.138 + 0.012 1445 5.43 +0.17 4.77+0.19 Kim et al. (2021) 
KMT-2017-BLG-2820 0.288 + 0.015 1.096 + 0.079 14.31 7.05 +0.44 5.94 +0.37 Ryu et al. (2021) 
OGLE-2012-BLG-1323 0.155 + 0.005 5.03 + 0.07 14.09 11.9+0.5 — 2.37 X 0.10 Mróz et al. (2019b) 
OGLE-2016-BLG-1540 0.320 + 0.003 1.65 +0.01 1351 151408 9.2 + 0.5 Mróz et al. (2018) 
OGLE-2019-BLG-0551 0.381 + 0.017 4.49 + 0.15 12.61 19.5 1.6 4.35 +0.34 Mróz et al. (2020c) 
MOA-9y-1944" 1.594 + 0.136 0.00928 + 0.00032 20.14 0.43 +0.10 46.1 + 10.5 K23 
OGLE-2017-BLG-0560^ 0.905 + 0.005 0.901 + 0.005 12.47 349 + 1.5 38.7 + 1.6 Mróz et al. (2019b) 


“Likely Brown dwarf lens. 


of the terrestrial mass population which gives rise to 
events such as OGLE-2016-BLG-1928 which has 0g = 
0.842 + 0.064 (Mróz et al. 2020b). These values are 
significantly smaller than the lower edge of 0g ~ 4.35 pas 
as reported in Gould et al. (2022). This is likely a result 
of selection bias given that Gould et al. (2022) focused 
on the sample with super-giant sources, see Figure 2. 

We compare the parameters of these events to 
seven known FFP candidates with 0g measurements 
in Table 1. The sources of all known FFP candi- 
dates except OGLE-2016-BLG-1928 are red clump 
giants (RCGs) or red super-giants which have large 
9,. = 5.4, 7.1,11.9, 15.1, 19.5 and 34.9 was. The magnifi- 
cation tend to be suppressed by large 0. with small 
0g, i.e., large p as Apgmax = wV1-4/p? (p > 1) 
(Maeder 1973; Agol 2003; Riffeser et al. 2006). For 
example, in case of the terrestrial mass lens with 
0g ~ lpas, the maximum magnification will be only 
Arg max = 1.066, 1.039, 1.014, 1.009, 1.005 and 1.002 for 
the above values of 0,, respectively. Note that the 
source of the terrestrial FFP candidate event, OGLE- 
2016-BLG-19288 is a sub-giant with 0, = 2.37 uas. It 
is important to search for short FSPL with sub-giants 
and dwarf sources to find low mass FFP. There is no 
FSPL event with a red super-giant source in our sample 
because these are saturated in MOA image data. 

K23 provides the detection efficiency for our FSPL 
event sample to compare to those in Gould et al. (2022) 
who tried to constrain the FFP MF from only g dis- 
tribution for their FSPL event sample. However, we 
perform a likelihood analysis using not only Og of the 
FSPL events, but also include tg of PSPL events which 
are more informative than Og distribution alone. See 
more details in Section 4. 


4. LIKELIHOOD ANALYSIS OF MASS FUNCTION 


In the final sample of K23, there are 10 (12) short 
timescale events with tg « 1 day after applying CR2 
(CRI). Figure 3 shows the tg distribution of the CR2 
sample. The distribution is roughly symmetric in log tp, 
with a tail at tg < 0.5. This confirmed the existence of 
such short timescale events with tg « 0.5 day as re- 
ported by Mróz et al. (2017). In this section, we per- 
form a likelihood analysis on each of the 3554 (CR1) 
and 3535 (CR2) events using a Galaxy model to con- 
strain the mass function of lens objects. 

We define the likelihood, £, in Section 4.1. In Sections 
4.2 and 4.3, we determine the mass function without and 
with a planetary mass population, respectively, by mini- 
mizing x? = —21n L. Although the absolute value of x? 
is not meaningful due to its dependence on an arbitrary 
normalization associated with our likelihood calculation, 
the fitting procedure is still statistically valid as the rela- 
tive likelihood between two models, represented by Ax?, 
is independent of the normalization. 

Note that results of the likelihood analysis for sample 
CR1 and CR2 are very similar. In the following sections, 
we show only the results for CR2 as our final results 
except in the tables. 


4.1. Likelihood 


Although our sample contains more than 3500 events, 
the mass function of planetary-mass objects is largely 
determined by the events with tg « 1 day, which ac- 
count for about 0.3% of them. We separately define two 
likelihoods; Lshort for short timescale events with the 
best-fit tg < 1 day, and Liong for events with the best- 
fit tg > 1 day. In our likelihood analysis, we use the 
combined likelihood of £ = Lshort Llong- 
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Figure 2. Extinction free CMD of gb3-7-6. The orange 
curve is the isochrone matched to this subfield. The cyan 
square is the RCG centroid. The red circles with error bars 
are sources of the 13 FSPL events in this work. The blue 
filled circles indicates the 2 FFP candidates in this work. 
The black open and filled circles are FSPL events and FFP 
events from Gould et al. (2022), respectively. The purple 
triangle indicates the source of terrestrial FFP, OGLE-2016- 
BLG-1928S (Mróz et al. 2020b). 


For Liong, we simply use the best-fit fg values pro- 
vided by K23, which is similar to the approach by pre- 
vious studies (Sumi et al. 2011; Mróz et al. 2017). This 
is because (i) the relatively smaller errors of tp, (ii) the 
effect of individual tg error is statistically marginalized 
by the large number of events, (iii) the limited sensitiv- 
ity to Og, and (iv) the minimal impact on our primary 
goal of measuring the mass function of planetary mass 
objects. 

On the other hand, the situation is the opposite for the 
short events, Lshort, ie., (i) the tg errors are relatively 
large due to their shorter magnification period, although 
smaller than the event selection threshold (see Table 2 
of K23), (ii) the number of events is very limited (12 for 
CR1 and 10 for CR2), and tg < 1 day is only sparsely 
covered in Figure 3. Thus, these may not be sufficient 
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Figure 3. The observed timescale tg distribution passing 
criteria CR2 from the 9 year MOA-II survey.. The red line 
indicates the best fit single lens model for all population. The 
blue dotted line represents the known populations of stars, 
brown dwarfs, and stellar remnants, and the green dashed 
line represents the planetary mass population. 


to statistically marginalize the effect of tg errors of in- 
dividual events in the likelihood analysis, (iii) because 
p = 0,/0g are expected to be larger than those of longer 
timescale events, one may get beneficial constraints on 
0g even if Og are not well determined, and (iv) they play 
a crucial role in determining the mass function of plane- 
tary mass objects. Therefore, we use the joint posterior 
probability distribution of (tg, 0g) for each event derived 
by K23 using the Markov Chain Monte Carlo (MCMC) 
method for Lsnort, to take into account both tg and 0g 
information and their uncertainties. 

We first describe the simpler one, Liong, in Section 
4.1.1, then describe Lshort in Section 4.1.2. 


4.1.1. Likelihood for events with tg > 1 day 
We define the likelihood for events with tg > 1 day by 


Niong 


Liong x II C (tg i; D), (4) 


i=1 


where i runs over all the Niong events that have the best- 
fit tg > 1 day in our sample (Mong = 3542 for CR1 and 
Mong = 3525 for CR2), and tg, is the best-fit tg value 
for ith event given by K23. 

The function G(tg;T) is the model’s detectable event 
rate as a function of tg with given model event rate I, 
combined for the 20 survey fields, given by 


G(tg; T) = 3 w; gj (te; T). (5) 
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Here, 7 takes field index values gb1 to gb21, except for 
gb6. See Table 1 of K23 for the location and properties 
of each field. The weight w; for the jth field is given by 


wj = » nRo,kfLF,k» (6) 


kej 


where k indicates a 1024 pixel x 1024 pixel subframe 
in the jth field (k = 1,2,...,80), nrc, is the number 
density of RCGs in the kth subfield, fry; is the fraction 
of stars with magnitude I < 21.4 mag in the kth subfield, 
and wj; is thus proportional to the expected event rate 
in the jth field. To calculate fLF,k, we used a combined 
luminosity function that uses the OGLE-III photometry 
map (Szymański 2011) for bright stars and the Hubble 
Space Telescope data by (Holtzman et al. 1998) for faint 
stars. 

The function g; is the model's detectable event rate 
as a function of tg for field j as given by 


gj (te; V5) = €j(te; D5) D; (te), (7) 


where é(tg;D) is the integrated detection efficiency of 
the survey as a function of tg. K23 demonstrated that 
when finite source effects are important, the detection 
efficiency, e(tg, 0g) is a function of two variables, tg and 
0g. Therefore, we must integrate over 0g to obtain the 
integrated detection efficiency, €(tg;I), which now de- 
pends upon the event rate and the mass function of the 
lens objects. This gives 


€j (tg; I) =| ej(tn, 0g) Dj (On|te) don, (8) 


where ej(tg,0g) is the detection efficiency for events 
with tg and Op for ith field. We use the detection ef- 
ficiency ej(tg, 0g) estimated by the image level simula- 
tions in K23 for the 20 fields of the MOA-II 9-yr survey. 

We consider the model event rate as a function of tg 
and (tg, 0g), l'(tg) and T (tg, 0g), respectively, as a nor- 
malized function so that its integration gives one, i.e., 
these are probability density functions of tg and (tg, 0g), 
respectively. D(0g|tg) = I(tg,0g)/T(tg) is the prob- 
ability density of events with 0g given tg. Thus, the 
calculation of €(tg;T) in Eq. (8) has to be done for ev- 
ery proposed MF during the fitting procedure because 
I'(0g|tg) depends on the MF. 

The function L';(tz, 0g) for jth field can be separated 
from the MF (Han & Gould 1996), 


D;(tg,0g) = J ste osa yegr raa, 
(9) 


where 7;(tz, 9m) is the event rate for lenses with mass 
1 Mo and (M) is the present-day MF (expressed as 
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Figure 4. Integrated detection efficiencies, (tp; DI), as a 
function of the timescale tz down to the source magnitude 
of Is < 21.4 mag for the criteria CR2. Red, black, green and 
blue lines indicate the efficiencies of fields with the highest, 
high, medium and low cadence, respectively. 


dN/dM). Although substituting Eqs. (8) and (9) makes 
the calculation of g;(tg;DI;) in Eq. (7) a double integral 
over M and 0g, K23 showed that the integration over 0g 
is largely avoidable during a fitting procedure by switch- 
ing the order of the integrals and calculating the integral 
over Og before the fitting. 

We calculate 7;(tz, 9m) for each field using the den- 
sity and velocity distribution of stars from the latest 
parametric Galactic model toward the Galactic Bulge 
based on Gaia and microlensing data (Koshimoto et al. 
2021a). 

Figure 4 shows the integrated detection efficiencies 
€(tg; D) for the event rate calculated with the best fit 
MF model with the criteria CR2. The curve for CR1 is 
similar. 


4.1.2. Likelihood for short timescale (tg < 1 day) events 


We follow Hogg et al. (2010) to utilize the posterior 
probability distribution of each event from MCMC to 
calculate the likelihood for the short timescale events, 
Lshort- Given the output MCMC samples of posterior 
distributions for individual events by K23, the likelihood 
is presented by 


Nshort Ki 


L short e X 


` G(log teik, log 8g, ik; D) 
1—1 k=l 


po(log teik, log kik) 


(10) 


where i runs over all the Nshort events that have the 
best-fit tg < 1 day (Nehort = 12 for CR1 and Nshort = 
10 for CR2), k runs over all the K; samples in the 
MCMC sample of posterior distribution for ith event, 
and po(log tg, log 0g) is the prior distribution multiplied 
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in the posterior. The model’s detectable event rate as a 
function of (log tg, log 0g) is given by 


G(log ty, log 8g; T) = S ^w; 9; (log tn, log @n;T';) (11) 
j 


with 


gj(log tg, log 0g; T4) = €; (log tp, log 0g) T; (log tg, log Og), 
(12) 


where we represented it as a function of (log tg, log 0g) 
rather than (tg, 0g) because K23 provide the posterior 
distributions with the uniform prior in (log tg, log 6g), 
i.e., po(log tg, log 0g) = const.. 

Eq. (10) calculates the likelihood by summing the ra- 
tio of G (log tg, log 0g; T) to po(log tz, log 0g) to replace 
the used prior (i.e., po) with the new prior (i.e., G) over 
the MCMC samples. This method, which uses all the 
MCMC samples, allows Lshort to account for the uncer- 
tainty of the parameters, unlike Llong given in Eq. (4). 

Despite the significant computational cost of Eq. (10) 
associated with performing a summation over K; (typi- 
cally ~ 5x 10°) samples for each proposed mass function 
during the fitting process, we addressed this by imple- 
menting a binning strategy for the MCMC sample using 
grids of (log tg, logg) with a size of (0.05 dex x 0.05 
dex), which significantly increased the computational ef- 
ficiency. 


4.2. Mass function of known population 


Firstly, we perform the likelihood analysis without the 
short events with tg « 1 day using the Galactic model 
with the MF of known population, i.e., stellar rem- 
nants (black holes (BH), neutron stars (NS) and white 
dwarfs(WD)), main sequence stars (MS) and brown 
dwarfs (BD). We use a broken power-law MF given by 


M-* (Mi < M/Me < 120) 


n^ c 4 M-? (0.08 < M/Mo < Mı) 
M-?* (3x107* < M/Mo < 0.08). 
(13) 
We adopt the values of parameters a; = 1.32 and 


a2 = 0.13, a3 = —0.82 and M, = 0.86 from the E+Ex 
model of Koshimoto et al. (2021a) by default unless 
specified as fitting parameters in the following three 
models. The minimum mass 3 x 1074 Mọ is taken to be 
smaller than the theoretical minimum mass of the gas 
cloud, ~Jupiter-mass, that collapses to form a brown 
dwarf (Boss et al. 2003). During our fitting procedure, a 
proposed initial mass function (IMF) is converted into a 
present-day mass function following the procedure used 


by Koshimoto et al. (2021a) that combines their stel- 
lar age distribution and the initial-final mass relation by 
Lam et al. (2020) to evolve stars into stellar remnants. 

We consider three models here: BD1, BD2, and BD3. 
In BD1, we fit only ag as a fitting parameter, while fix- 
ing o4, a», and Mj. Similarly, in BD2, we fit a3 and 
0o», and in BD3, we fit o1, a2, ag, and Mi. To per- 
form the fitting, we use the Markov Chain Monte Carlo 
(MCMC) method (Metropolis et al. 1953), and assign 
uniform distributions as priors for all the parameters. 

'The best fit models BD1, BD2 and BD3 are almost 
indistinguishable from the blue dotted line in Figure 3. 
One can see that the models fit the data with tg > 1 
day very well. The best fit parameters and x? values 
are listed in Table 2. There is no significant difference 
in the resultant parameters between different selection 
criteria or among the BD1, BD2, and BD3 models. 

All of the parameters are consistent with those of 
Koshimoto et al. (2021a) within lc. This indicates 
that our dataset confirmed the Galactic model and MF 
of known objects by Koshimoto et al. (2021a). This 
also indicates that our dataset is consistent with the 
OGLE-IV tg distribution for tg > 1 day (Mróz et al. 
2017, 2019) that is fitted by Koshimoto et al. (2021a). 

In the following analysis, we fit only o3 and fix all 
other parameters for the known populations. Note, in 
Koshimoto et al. (2021a), the Galactic model and MF 
are constrained to satisfy the microlensing tg distribu- 
tion, stellar number counts and the Galactic Bulge mass 
from other observations, simultaneously. In principle, 
the MF should not be changed alone because it is related 
to other parameters of the Galactic model. However, the 
contribution of objects with M/Mọo < 0.08 are negligible 
in stellar number counts and as a fraction of the Galac- 
tic Bulge mass. Thus, we assume that a model with a 
different slope at lower masses with M/Mo < 0.08 is 
still valid. 


4.3. Mass function of planetary mass population 


If the candidates with tg < 0.5 day are really due to 
microlensing, they can not be explained by known popu- 
lations, i.e., stellar remnants, MS or BD. To explain the 
tail for short values of tg, we defined a new model “PL” 
which introduces a planetary mass population by the 
following power law in addition to known populations 
(Eq. 13), 


Mnorm 


dN, -z M 
dlog M 


E 
) (Mais € M/Mo < 0.02). 

(14) 
Here Z is a normalization factor and Myorm is a ref- 
erence mass whose inclusion allows Z to have a unit 
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Table 2. Best fit parameters of the mass function for known population. 


model BD1 BD2 BD3 Koshimoto+21la* 
CR1 CR2 CR1 CR2 CRI CR2 

Mi (0.86) (0.86) (0.86) (0.86) 0:970 04 0.99 089 0.86* 0 00 

oa (1.32) (1.32) (1.32) (1.82). dX38t037. — q:34t02 1.327015 

az (0.13) (013) | 0.200007 =o. 20F8- 62  g23t795 — 0.2470-07 0:18*012 

a3 —0.60*0 09 0.627099 0.74913  —Qv6*0 6 0.7619 — —0.79*0:32 —0.82*0 21 

x? 35919.4 35722.6 35918.2 35721.5 35918.0 35721.3 


“Results of fitting to various bulge data including the OGLE-IV tg distribution of tg > 1 day (Mróz et al. 2017, 
2019). The representative values are shifted to the ones for the E+Ex model from their original ones for the 


G+Gx model. 


NoTE—Some of the upper errors of Mı is negative because the best fit value is outside of the 68% range. This is 


because M, is restricted to be less than 1 Mo. 


of (dex)^!. Although Myorm can be an arbitrary zero 
point, we found that the uncertainty in Z is minimized 
when we adopt Mnorm = 8 Me which is recognized as a 
pivot point. 

In the model PL, we use o3, o4 and Z, as fitting pa- 
rameters and fix parameters o, = 1.32, a» = 0.13 and 
Mı = 0.86 (Koshimoto et al. 2021a). We assign uni- 
form distributions as priors for a3, a4 and log Z in our 
MCMC run. We found that the fitting result does not 
depend on Mmin at all when Mmin < 3x 1077 Mo, which 
indicates our data sensitivity is down to ~ 3 x 107" Mo. 
Thus, we decided to use Mmin = 1077. 

'The red solid line in Figure 3 represents the best fit 
model for all populations with the CR2 sample. This 
figure indicates that the model represents the observed 
tg distribution well. Note that although the observed 
tg distribution shown in black in Figure 3 does not in- 
clude error bars along the tg axis, the best-fit line is 
derived from our likelihood analysis that takes into ac- 
count the tg errors as well as the Og constraints for the 
short events with tg « 1 day. Figure 5 shows the pos- 
terior distributions of the parameters of PL model. The 
best fit parameters and x? are listed in Table 3. 

The best fit power index for BD is o3 = —0.58* 0-12 
which is consistent with the model without the planetary 
mass population. 

The best fit MF of the planetary mass popu- 
lations with the normalization Z relative to stars 
(MS+BD+WD) (integrated IMF over 3 x 107^ < 
M/Mo < 8) can be expressed as 


N, 2.187952 M —O4 
d 4 —1.40 ( ) (15) 


dlog M ~ dex x star 8 Ma 


where a4 = 0.96*0 51. Figure 6 shows the IMF of the 
best fit PL model. This o4 is consistent with the corre- 


sponding power law index of 0.9 S p S 1.2 reported by 
Gould et al. (2022). 

'This can be translated to the normalization per stellar 
mass of stars, Z Mo, as, 


N. 5.481118 / aN 
dN4 —3.50 ( ) (16) 


dlogM | dex x Mo \8 Mẹ 


This implies that the number of FFPs per stars is 
f= PU star! over the mass range 109 < M/Ms < 
0.02. Note taht this value is vary depending on the 
minimum mass. The total mass of FFPs per star is 
m = 80113 Mo (0.25+9:23 Mj) star-!. This is less depen- 
dent from the minimum mass. The total mass of FFPs 
per Mo is mMe = 2024198 Ma (0.64*0.1? M;) Mg. This 
is more robust values less dependent on uncertainty in 
the abundances of the low mass objects for both FFP 
and BD. 

'The normalization, number and total mass of FFP 
relative to MS+BD (3 x 10^ < M/Mo < 1.1) are 
also shown in Table 3. These normalizations can be 


translated to Zwsipp = 0.531919 dex-!star-! and 


Zuo. gp = 244191 dex-1 M! with Mnorm = 38Mọ. 
These are almost same as Zygipp = 0.39 + 0.18 
dex-!star-! and Zms+BD = 1.964£0.98 dex-! M! with 
Myorm = 38M by Gould et al. (2022). 

Note that the lenses for these short events could be 
either FFP or planets with very wide separations of more 
than about ten astronomical units (AU) from their host 
stars, for which we cannot detect the host star in the 
light curves. 


5. DISCUSSION AND CONCLUSIONS 
We derived the MF of lens objects from the 9-year 
MOA-II survey towards the Galactic Bulge. The 3,535 
high quality single lens light curves used in our statistical 
analysis include 10 very short (tg < 1 day) events, and 
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Table 3. Best fit parameters of the mass function for the planetary 


mass population. 


3r 2 
x xBEST AX,<1 
S AX <4 
te AX <9 
T Ay®<16 
Ay*=16 
N o 
op 
© 
T -1 
1.5x10* L 
4 
z, 10 - 
5000 E 
0 
-1.5 -1 -05 0 1 2 3 
[01 3 Oy 
Figure 5. Posterior distributions of the parameters of the 


PL model for sample CR2. The vertical red dotted lines indi- 
cate the median and +1ø. The vertical orange line indicates 
the best fit. 


Mass [Mo] 
107 10? 10? 10! 10! 
F | T | | 
103 N —— All 
E == Planet 
102 L zeae Star+BD 
E —-— Bound planet 
10? F 
L E 
= E 
10° $ 
1071 $ 
10-2 E 
E | l | l 
10° 10? 10^ 109 
Mass (Mo) 
Figure 6. Initial mass function (IMF) of the best fit PL 


model for CR2. The red line indicates the best fit for all 
population. The blue dotted line and green dashed line show 
the IMFs for the stellar and brown dwarf population and 
for the planetary mass population, respectively. The shaded 
areas indicate lo error. The gray dashed-line and the shaded 
area indicate the best-fit and 1 o range of the bound planet 
MF by Suzuki et al. (2016) via microlens. 


CR1 CR2 Gould+22 
(Mnorm) (8 Mg) (8 Mo) (38 Mg) (38 Mg) 
Mi (0.86) (0.86) 
o1 (1.32) (1.32) 
az (0.13) (0.13) 
as —0.55*012 — —0.58*0 15 
as 0.90* 0:25 0.96* 9:27 fixed at 0.9 or 1.2 
+0.54 +0.52 +0.17 
Z 2.087 133 2.187 145 0.49 5.37 
ZMs+BD 227 noe 2.881099 0.5345 18 0.39 + 0.204? 
zMo 5.83: 126 BAgt EE. — 1.22025 
M, E 
Z IS UHD 10.638252 — 10.905299 2.44* 0 0) 1.96 + 0.98? 
a +20 +23 
f 1730 21125 
fus+B0° 1975 23775 
Maa +54 +59 
ee 457 30 53134 
© a 107 +117 
fus+Bp 89* 59 106258 
mP 89* 96 Mo 80* 72 Me 
mMS- BD? 98* 17 Mo 88* 91 Me 
m™oP 2297718Me 2027199 Me 
M 
musieo? 4571339Me 4041333 Mo 
x? 36273.0 36024.1 


“Number of planetary mass objects per BD--MS--WD (f), per MS+BD 
(fms+BD), per solar mass of BD+MS+WD (f™“©) or per solar mass of 
MS+BD (fi. 


gO) 
Ius+Bp 


) when MF down to 


vary depending on the minimum mass. 


0-8 Mo are integrated. These are 


P Total mass of planetary mass objects per BD+MS+WD (m), per MS+BD 
(mms+BD), per solar mass of BD+MS+WD (m9) or per solar mass of 


MS+BD (my pp) when MF down to 1079 Mo are integrated. 


NoTE—We adopt the model for CR2 as the final result. 


13 events with strong finite source effects that allow the 
determination of the angular Einstein radius, 0g. 

The cumulative 0g histogram for these 13 events re- 
veals an “Einstein gap" at 5 < Og/pas < 70 which is 
roughly consistent with the gap at 10 < Op/pas < 30 
found by the KMTNet group (Ryu et al. 2021; Gould et 
al. 2022). This gap indicates that there is a distinct 
planetary mass population separated from the known 
populations of brown dwarfs, stars and stellar remnants. 

We constructed the tg distribution of all selected sam- 
ples including both PSPL and FSPL. We calculated 
the integrated detection efficiency &(tg; D) of the survey 
by integrating the two dimensional detection efficiency, 
€(tg, 0g), measured from image level simulations that in- 
cluded the FS effect, and convolving this with the event 
rate I'(tg, 0g) given by a Galactic model and MF. We 
found that the tg distribution has an excess at short 
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tg values which can not be explained by known popula- 
tions. 

We then adopted the power law MF for the plane- 
tary mass population. We found that these short events 
can be well modeled by dN4/dlog M = (2.18193) x 
(M/8Mg)-?* dex star! with ag = 0.96103! at 
1077 < M/Mg < 0.02 (or 0.033 < M/Ma < 6660). 

This can also be expressed by the MF per stellar 
mass as, dN4/dlogM = 5.48*115 x (M/8Mg)-?* 
dex-! M5'. We showed the number of FFP or distant 
planets is f — 21733 per stars. 

It is well known that planet-planet scattering during 
the planet formation process is likely to produce a pop- 
ulation of unbound or wide orbit planetary mass ob- 
jects (Rasio & Ford 1996; Weidenschilling & Marzari 
1996; Lin & Ida 1997) The probability of planet scatter- 
ing likely increases with declining mass because planets 
usually require more massive planets to scatter. So, we 
expect the power law index of MF of bound planets ap 
is smaller than that of a4 for unbound or large orbit 
planets, i.e., a4 > ap. 

One can compare our FFP result to the MF of known 
bound planets. At present, microlensing surveys have 
only measured the mass ratio function, rather than the 
mass function, of the bound planets. Currently, the 
most sensitive study of the bound planet mass ratio 
function Suzuki et al. (2016) found that the mass ra- 
tio function can be well explained by the broken power 
law with a, = 0.93 + 0.13 for q > qbr = 1.7 x 1074, 
ap = —0.6*03 for q < gr = 1.7 x 1074. While the 
Suzuki et al. (2016) data could establish the existence 
of the power-law break with reasonably high confidence 
(a Bayes factor of 21), there was a large, correlated un- 
certainty in the mass ratio of the break and slope of the 
mass ratio function below the break. So, we chose to fix 
the mass ratio of break at qbr = 1.7 x 1074 in order to 
estimate the power law below the break. 

More recently, several papers have attempted to im- 
prove upon this estimate by including a heterogeneous 
set of lower mass ratio planets found by a number of 
groups without a calculation of the detection efficiency. 
These efforts included attempts to estimate the effect of 
a “publication bias” that might cause planets deemed to 
be of greater interest to be published much more quickly, 
leading to biased, inhomogeneous sample of planets. 
This “publication bias” is caused by the decision to pub- 
lish some planet discoveries at a higher priority than oth- 
ers. With such an analysis Udalski et al. (2018) reported 
ap = —1.05*0$8 with their sample and ay = —0.7319°34 
when combined with the Suzuki et al. (2016) result for 
q<1x 10-4 < qr. A similar analysis by Jung et al. 
(2019), attempted a new measurement of the location of 


the break and found ap = —4.5 for q < qpr = 0.55 x 1074 
which is consistent with the Suzuki et al. (2016) result 
when qbr is not fixed. However, a more recent paper 
(Zang et al. 2022) by many of the same authors, re- 
ported a number of planetary microlensing events that 
were missed by the analyses described in Udalski et al. 
(2018) and Jung et al. (2019). This casts some doubt 
on the validity of some of the assumptions in these pa- 
pers. This later paper also suggests that planets with 
mass ratios of q < qbr = 1.7 x 107+ may be more com- 
mon than previously thought, although a more definitive 
claim awaits a detection efficiency calculation. Also, the 
Suzuki et al. (2016) analysis does not imply that there 
is a peak in the mass ratio. Instead it concludes that 
the slope does not rise as steeply toward low mass ratios 
as is does for q > 1.7 x 1074. 

The broken power-law model of Suzuki et al. (2016) 
is consistent with the hypothesis that these unbound or 
wide orbit planetary mass objects are the result of scat- 
tering from bound planetary systems. It is the lower 
mass planets that are preferentially removed by planet- 
planet scattering interactions, so the initial planetary 
mass function may have been closer to a single power-law 
with ap ~ 0.9, but planet-planet scattering has likely de- 
pleted the numbers of low-mass planets at separations 
beyond the snow line where microlensing is most sen- 
sitive. Thus, planet-planet scattering may be responsi- 
ble for the mass ratio function “break” observed in the 
Suzuki et al. (2016) sample 

'This idea that planet-planet scattering is responsible 
for a FFP mass function slope that is steeper than the 
slope of the mass ratio function for low-mass bound 
planets is also consistent with the single power-law mod- 
els that were found in smaller data sets (Sumi et al. 
2010). The best fit single power-law model for the 
Suzuki et al. (2016) sample gives dNpouna/dlogg = 
0.068* 6016 dex-?star-! x (q/0.001)-9* with os, = 0.58+ 
0.08 for 3x 1079 < q < 3x 107?, but the broken power- 
law is a significantly better fit to the Suzuki et al. (2016) 
data. Note, this single power law model with ap = 0.58 
satisfies a4 > ap, for our value of o4 = 0.96* 027, imply- 
ing that unbound (or very wide orbit) planets increase 
more rapidly than bound planets at low masses. Thus 
our main conclusion discussed bellow with the broken 
power law model, which the lower mass planets are in- 
creasingly scattered, is not specific to the Suzuki et al. 
(2016) broken power-law model. 

As a comparison, we transformed the bound planet’s 
mass “ratio” function of Suzuki et al. (2016) to a mass 
function by using the estimated average mass of their 
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hosts of ~ 0.56M as shown? in Figure 6. We estimate 
the abundance of the wide-orbit bound planets to be 
fwide = LITE planets star™t in the mass range 1076 < 
M/Mo < 0.02 (0.33 < M/Mẹ < 6660) and separation 
range 0.3 < s < 5, which corresponds to a semi-major 
axis of roughly 0.7 < a/au < 12. This indicates that the 
abundance of FFP, f = apes planets star-!, is 19773 
times more than wide-orbit bound planets in this mass 
range. 

This is because the number of wide-orbit bound plan- 
ets decreases at lower masses than the break at Mpreak © 
1.0 x 1074Mọo, while the number of high-mass bound 
planets is larger than that for FFP. Again, this is con- 
sistent with the hypothesis that the low-mass planets are 
more likely to be scattered. Note that there is still large 
uncertainty in the MF at low masses for both bound and 
unbound planets. It is very important to constrain the 
these MFs at low masses. 

We can also compare our number for the FFP abun- 
dance with the abundance of the bound planets with 
short period orbits of P — 0.5 — 256 days and planetary 
radii of Rp = 0.5 — 4Rẹ found by Kepler. Hsu, Ford, 
Ragozzine & Ashby (2019) find frak = 3.519% for FGK 
dwarfs and Hsu, Ford & Terrien (2020) find fm = dre 
for M dwarfs. Because the typical spectral types of their 
samples are G2 (M = 1Mọ) and M2.5 (M = 0.4Mo), 
their typical semi-major axis are 0.012 < a/au < 0.79 
and 0.009 < a/au < 0.58, respectively. The fraction 
of FGK and M dwarfs relative to all population except 
BH and NS are 0.157 : 0.465 in our best fit MF. By 
weighting with these stellar type fractions, the abun- 
dance of the known close-orbit bound planets is about 
felase. = 25005 per star. (This ignores the relatively 
small number of gas giant planets in short period orbits 
(Bryant et al. 2023)). 

The total abundance of the wide-orbit and known 
close-orbit bound planets is about fpouna = 3-670) 
per star. This indicates that the abundance of FFP, 
f= giu planets star-!, is E Mee times more than 
known bound planets in this mass range. 

We found the total mass of FFPs or distant planets 
per star is m = 80*722M&(0.25*0 22 My) star-! in this 
107 < M/Mọ < 0.02 (0.33 < M/Ma < 6660) mass 
range. This is comparable to the value of 91733 Ma 
star-! for wide-orbit bound planets with separations of 
0.3 < s < 5 in the same mass range. It is not straight 
forward to estimate the total mass of inner planet found 


? The lo range indicated by the gray shaded area in Figure 6 
does not match the one provided in Suzuki et al. (2016). This was 
due to an error in the Suzuki et al. (2016) figure, but there is no 
error in the other results in that paper. 


by Kepler because only a small, and somewhat biased, 
sample of Kepler planets have mass measurements. The 
total masses of FFP and bound planets are less depen- 
dent on the uncertainty of the number of low mass plan- 
ets than the total numbers of FFP and bound planets 
are. 

These comparisons indicate that 1919 times more 
planets than the ones currently in wide orbits have been 
ejected to unbound or very wide orbits. These com- 
parisons also suggest that the total mass of scattered 
planets is of the same order as those remaining bound 
in wide orbits (beyond the snow line) in their planetary 
systems. The low mass bound planets in wide orbits are 
much less abundant than those orbiting closer to their 
host stars. This may be explained by that planets in 
wide orbits are more easily ejected than those in close 
orbit. 

'The power-law index of the IMF of planets formed 
in wide orbits in protoplanetary disks is likely to be 
a4 ~ 0.9 with an abundance of 29 planets star-! 
or 171189 Me (0.54*029 M3) star—!. 

Another, rather speculative, possibility is that most 
of the low-mass objects found by microlensing are pri- 
mordial black holes (PBH). Hashino et al. (2022) pre- 
dicted PBH generated at a first order electroweak phase 
transition have masses of about 107? Mo. They found 
that depending on parameters of the phase transition 
a sufficient number of PBH can be produced to be ob- 
served by current and future microlensing surveys. The 
mass of such PBH is a function of the time of their 
generation, i.e., the electroweak phase transition, and is 
expected to be a delta-function distribution. To differ- 
entiate PBH from FFP, we need to measure the shape 
of the MF accurately. This can be done by the cur- 
rent (MOA, OGLE, KMTNet) surveys, the near future 
(PRIME) ground telescope and the Roman Space tele- 
scope. 

For the first time, we have determined the detection 
efficiency as a function of both the Einstein radius cross- 
ing time and the angular Einstein radius, because finite 
source effects have a large influence on the detectability 
of microlensing events due to low-mass planets. This 
method is necessary for reliable results for low-mass 
FFPs, and it should be very useful for the analysis of 
these future surveys which will detect many short events. 

A precise measurement of the free floating planet mass 
function will require a microlensing survey that can ob- 
tain precise photometry of main sequence stars with 
relatively low magnification, because the small angular 
Einstein radii, 0g, of low-mass planetary lenses prevent 
high magnification. The exoplanet microlensing survey 
of the Roman Space Telescope is such a survey, and it 
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should provide the definitive measurement of the free 
floating planet mass function. 
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